Modified shelfwave_set_OBC_data to work properly with grid rotation of 90, 180
or 270 degrees.  With these changes the ESMG shelfwave test case is now giving
bitwise identical answers under grid rotation.  This commit also corrects a
recently added problem with potentially unallocated arrays being used to test
whether a pointer is set to the unrotated tracer reservoirs in cases where there
are no tracers.  This commit does change (and fix) answers when grid rotation is
used, but answers are bitwise identical when there is no grid rotation, and
there are no changes to any public interfaces.

  Modified tidal_bay_set_OBC_data to work properly with grid rotation of 90, 180
or 270 degrees.  With these changes the ESMG Tidal_bay test case is now giving
bitwise identical answers under grid rotation.  This commit does change (and
fix) answers when grid rotation is used, but answers are bitwise identical
when there is no grid rotation, and there are no changes to any public
interfaces.

* *Implement subcycing for internal tides refraction/propagation

This PR allows to bring back the thermodynamics timestep higher
while keeping the raytracing under CFL by doing multiple steps
within the thermo step.

The only variable that needs to be accumulated in the subcycling
is the slope loss (residual of transmission and reflection on slopes)
Since it is already accumulated within the 2 directions of propagation,
it only needs to be initialized outside of the subcycling.

* unindent loop for easier git diff (will indent back later)

* fix bad placement of complete_group_pass and sync En

* rewrite subcycle in terms of a timestep and group_pass

* fix indexing bug, remove useless init/passing

* fixes from review

---------

Co-authored-by: Raphael Dussin <Raphael.Dussin@noaa.gov>

  This commit adds the new runtime parameter ENABLE_BUGS_BY_DEFAULT and use it
to set the defaults for a number of bug-fix parameters that have been added
within about the past year and have not yet had their defaults set to correct
the bugs. The default value (True) allows us to move the version of the code
forward without changing any existing answers by default.  Conversely, setting
this to false means that no bugs are used unless they are actively selected,
which is probably the right choice for new runs.  The defaults for groups of
bug-fix flags are periodically changed to correct the bug, at which point this
new parameter will no longer be used to set their default.  This commit also
uses ENABLE_BUGS_BY_DEFAULT to set the defaults for 16 existing runtime
parameters, including  DETERMINE_TEMP_CONVERGENCE_BUG, RHO_PGF_REF_BUG,
VISC_REM_BUG,  OBC_RESERVOIR_INIT_BUG, OBC_HOR_INDEXING_BUG, EXTERIOR_OBC_BUG,
VERTEX_SHEAR_VISCOSITY_BUG, VERTEX_SHEAR_OBC_BUG, EPBL_MLD_ITER_BUG,
MEKE_GM_SRC_ALT_SLOPE_BUG, HOR_DIFF_LIMIT_BUG, MIXING_COEFS_OBC_BUG,
FRICTWORK_BUG, LA_MISALIGNMENT_BUG, IDL_HURR_SCM_EDGE_TAPER_BUG and
CHANNEL_FLOW_OBC_TRANSPORT_BUG, which are scattered across 16 modules.  For 7 of
these parameters, the default will be reset to false soon.  By default all
solutions are bitwise identical, but there is a new runtime parameter that will
occur in all of the MOM_parameter_doc.all files.

  Corrected the previous PE-count and layout dependency of answers for cases
with open boundaries are not at the natural edges of the domain by storing the
previous barotropic velocities in ubt_prev and vbt_prev over a wider range of
loop indices, and by increasing the stencil size used in btstep_timeloop in such
cases.  Because the required loop range no longer matches the range over which
the new velocities are set in bt_loop_update_u() and bt_loop_update_v(), the
code setting these arrays has been moved out of these routines and into the main
body of btstep_timeloop().  This commit will change answers (and make them
invariant to the PE count and layout) in some cases, such as the ESMG-configs
Kelvin_wave/rotate45_BT and Kelvin_wave/rotate_BT test cases, where there are
southern open boundaries appearing in the northern halo regions on any PEs, for
example.  In all the Kelvin_wave test cases, the answers now agree with the
previous single-PE answers.  Most regional configurations with OBCs only apply
them at the edges of the domain and hence will not be impacted, and no answers
are changed in cases without open boundary conditions.

  This commit adds the new runtime parameters BT_WIDE_HALO_MIN_STENCIL and
DEBUG_BT_WIDE_HALOS to control the stencil width used with wide halos and to
specify whether the checksums for the entire active halo of values are written
when DEBUG_BT is enabled, or just the symmetric computational domain.  Setting
`BT_WIDE_HALO_MIN_STENCIL = 2` and `DEBUG_BT_WIDE_HALOS = False` provides for
cleaner checksums in some cases because they are invariant to the use of wide
halos, the width of those wide halos, or the relative positions of interior-
domain OBCs relative to PE boundaries.  This commit also changes a number of
velocity point checksums in btstep_timeloop to do symmetric output.  All model
solutions are bitwise identical regardless of how these new parameters are set,
but they can alter some debugging output.

Inconsistencies in jobs within a workflow were leading to differences in
compiler versions across jobs.  This was breaking GitHub Actions
workflows.  We want to continue using `*-latest` runners, so we have
merged the FMS stage into each of the MOM6 builds.

There are additional changes which improve the overall workflow/job
decomposition.

* FMS builds included in MOM6 stage

* Artifacts are now packed into tarballs, which preserve most of the
  file state (mtime, file mode, etc).

  * No longer any need to chmod the executables.

  * `*.gnco` timestamp warnings are removed.

  * Overall simplification of generated executable lists; we now
    assume that build/* only contains executables.

* build.timing and build.unit are now correctly tagged as .NOTPARALLEL

* Macros for defining the target URL are now environment variables,
  stored in the workspace `env`.

  Fixed spatial indexing bugs with Kelvin_set_OBC_data when the new runtime
parameter KELVIN_SET_OBC_INDEXING_BUGS is set to false.  All answers are
mathematically equivalent even when KELVIN_SET_OBC_INDEXING_BUGS is true, but
only because in the existing cases using Kelvin_set_OBC_data use an f-plane with
a uniform depth, but there are changes at roundoff with some compilers for some
of the ESMG/Kelvin_wave test cases.  Also added comments noting where there are
missing temporal variability factors.  By default, all answers are bitwise
identical but there is a new runtime parameter.

  Modified Kelvin_set_OBC_data to work properly with grid rotation of 90, 180 or
270 degrees.  With these changes the ESMG Kelvin_wave test cases are now giving
(mostly) bitwise identical answers under grid rotation.  This commit does change
(and fix) answers when grid rotation is used.  There are no changes to any
public interfaces.

update MOM6 to its main repo. 20250527 commit.

  Do not set non-segment 3d OBC arrays in the unrotated OBC type when grid
rotation is in use, and added a fatal error if there is an attempt to use grid
rotation with tracer reservoirs when OBC_RESERVOIR_INIT_BUG is true.  Also
eliminated duplicative statements nullifying these pointers after they have been
deallocated.  Also converted the non-segment 3d OBC arrays from pointers to
allocatables, as there is no longer any reason for them to be pointers.  All
answers are bitwise identical unless undergoing a grid rotation test, and the
grid rotation tests work properly for all cases with the relevant bugs disabled.

  Added the new on_face element to the OBC_segment_data_type to indicate where a
field is staggered, and the used this to simplify the code and the logical tests
in initialize_segment_data and update_OBC_segment_data.  The new function
field_is_on_face() is used to set the on_face elements.  Also, duplicated code
that determines whether data requirements have been satisfied has been moved
into a single block of code that used regardless of whether the data is
specified or provided from a file.  There is also the new routine field_is_tidal
that can be used to determine the appropriate vertical extent of a field within
the model.  With these changes, there has been a substantial consolidation of
the blocks of code that allocates and initializes the various input buffers.

  In addition, some of the duplicated code setting the tracer reservoir values
for temperatures, salinities and other "obgc" tracers has been consolidated into
a single block of code.  The code setting the SSH at boundaries was also
consolidated rather than having duplicate blocks with and without the ramping
enabled.  Several redundant logical tests or unneeded arrays were also
eliminated.

  With these changes, MOM_open_boundary.F90 is shorter by about 190 lines and
should be easier to read and maintain.  There are no changes to publicly visible
interfaces, and all answers are bitwise identical.

  Fix a bug in determining where to apply OBC-related masking when interpolating
v-velocities to u-points to calculate the surface friction velocity under ice
shelves.  It would usually only arise at eastward OBCs under ice-shelves, but it
could also occur for westward OBCs that are within a few points of the western
side of the domain.  This bug breaks rotational symmetry and reproducibility
across processor counts and layout when there are any open boundary conditions
under ice shelves.  However, there are no reports of the detection of any such
problems, and I am unaware of any cases using open boundary conditions under ice
shelves, so it seems very likely that there are no existing MOM6 configurations
that would trigger this specific bug.  Given all these considerations and the
resulting expectation that it is exceptionally unlikely that fixing this bug
would impact any existing work, rather than introducing a bug-fix runtime
parameter to retain this bug, this commit simply fixes it with the obvious
single-character correction.

  Set the sign of OBC%segnum_u and OBC%segnum_v to indicate the direction of
open boundary currents at the u- and v- velocity faces, with positive values at
logically eastern and northern boundary condition faces and negative values at
western and souther boundary condition faces.  This new information is used to
set open boundary condition properties in 10 files, while in 5 of these files
the absolute value is used to recover the previous answers.  Also added the new
elements u_OBCs_on_PE and v_OBCs_on_PE to the ocean_OBC_type.  All answers are
bitwise identical, but there are two new elements and an important change to two
elements of the publicly visible transparent ocean_OBC_type.

  Add 16 integer elements to the ocean_OBC_type to indicate the ranges of
indices over which various directions of open boundary conditions are found in
the data and computational domain of the local PE.  Also added 4 logical
elements that indicate whether each of the 4 orientations of OBCs are fond on
the local PE.All answers are bitwise identical, but there are 16 new integer
elements and 4 new logical elements in a transparent public type.

  Refactored the OBC-related code in 5 routines (set_viscous_BBL(),
vertvisc_coef(), find_coupling_coef(), calc_isoneutral_slopes() and
ALE_remap_set_h_vel_OBC()) in 4 modules to use the restricted OBC loop ranges
for the various orientations of OBCs for greater efficiency, and in some cases
to avoid the use of multiple levels of nested types and looping over segments
that seem to be causing problems with our efforts to port MOM6 to GPUs.  The
cases where the loop over the OBC segments continue to be used are those where
the details of the segments matter (and not just their orientation) or where
they are being applied in ways where the order of the segments might matter.
All answers are bitwise identical, and there are no changes to any publicly
visible interfaces.

  Added the new parameter INTERPOLATED_SQG_STRUCTURE that can be set to false to
avoid interpolating properties to velocity points and then interpolating the
buoyancy frequencies and layer thicknesses back to tracer points when
calculating the SQG vertical structure.  This involves extensive additions to
calc_sqg_struct() to calculate the SQG vertical structure on an individual
column without interpolation.  The new code also cancels out common thicknesses
in the numerator and denominator and uses thickness_to_dz() to find the vertical
extents of the layers rather than using find_eta() and taking a difference
between the heights of pairs of interfaces, making it far less susceptible to
roundoff errors and able to work properly with vanishing layers.  The one change
when this new option is not set to false is that the diagnostics of N2_u and
N2_v are no longer posted when SQG_EXPO is not positive, although it is not
clear to me that these diagnostics should ever be posted from this routine.  The
default for INTERPOLATED_SQG_STRUCTURE is true, but I suspect that all of this
interpolation of N2 may be undesirable, especially as it might lead to
problematic behavior near coastlines or seamounts, and I will recommend that the
default be changed to false once this can be done gracefully without disrupting
any applications that are already using calc_sqg_struct().

  As a part of these changes, SQG_EXPO will only be logged when it would
actually be used, and the logic governing a number of deallocate calls in
VarMix_end() was simplified.  The Laplac3_const_u, Laplac3_const_v, KH_u_QG and
KH_v_QG elements of the VarMix_CS were converted from potentially static arrays
to allocatable arrays to reduce the MOM6 memory footprint when they are not
in use.  These should always have been allocatable, and there is no change
when MOM6 is not in static memory mode.

  By default all solutions are bitwise identical, but there is a new runtime
parameter that will appear in some MOM_parameter_doc files, while another unused
runtime parameter will disappear from others.

  Converted 28 arrays from potentially static using the ALLOCABLE_ and ALLOC_
macros into dynamic arrays in three modules.  This change reduces the static
memory footprint of MOM6.  These arrays are only ever used conditionally based
on parameter settings, so they should never have been declared using these
static memory macros.  All answers are bitwise identical and there are no
changes to any public interfaces.

* Add conditional halo pass for frazil inside make_frazil

Frazil may need an extra halo pass, but only when all of the following are true:
(1) frazil is calculated and applied in the halo (e.g., vertex_shear is true)
(2) the thermodynamic loop is called more than once before passing frazil to the sea-ice (dt_therm<dt_cpld)
(3) reclaim_frazil is true
(4) a previous make_frazil was called after the diabatic loop and before T&S were updated in halos.

This commit adds the halo pass inside make_frazil since it is where we can know all of these conditions at once.  We needed to add a flag to know if frazil was reset (which happens if dt_therm>=dt_cpld), in which case the halo pass would be unnecessary.  We needed an optional argument for make_frazil so we could flag to only do this update on the first call to make_frazil (which is needed because the previous call to make_frazil was after the diabatic loop, condition 4)

* Shorten line in MOM_diabatic_driver

* Move halo update on frazil out of make_frazil and into post_diabatic_halo_updates

---------

Co-authored-by: brandon.reichl <brandon.reichl@noaa.gov>

  Removed the option to call write_ocean_geometry_file() from within
MOM_initialize_fixed().  This option was only being exercised in one of the
three places where MOM_initialize_fixed was being called, and by doing so it
precluded the possibility of specifying an alternative name for this file.
Removing this option eliminates two arguments to MOM_initialize_fixed(), but it
requires the addition of a new call to write_ocean_geometry_file() directly from
shelf_main in ice_shelf_driver.F90.  All answers are bitwise identical, but
there is a simplification in the arguments to a publicly visible routine.

* Change MOM_generic_tracer into a stub

This PR removes the content from the functions and routines in MOM_generic_tracer to turn it into a stub.
In addition, the other modules in config_src/external/GFDL_ocean_BGC are removed since they are no longer called anywhere.

* Clean up module use statements

Remove unneeded use statements.

* Add fatal error if stub code is called

A fatal error has been added to register_MOM_generic_tracer which
will occur if the stub version of MOM_generic_tracer has been
compiled but use_generic_tracer is true.

* Fix missing '&' at continued line

---------

Co-authored-by: Theresa Cordero <Theresa.Cordero@gaea65.ncrc.gov>
Co-authored-by: Theresa Cordero <Theresa.Cordero@gaea63.ncrc.gov>

A new subroutine is added in MOM_interface_heights to calculate mass per
unit area [R Z ~> kg m-2] and bottom pressure [R L2 T-2 ~> Pa]. The
subroutine packages the original code for the diagnostics of these two
quantities in MOM_diagnostics and will be used for self-attraction and
loading.

A new option is added to calculate reference bottom pressure from the
initial state of the model run, based on tihckness, temperature and
salinity.

* The `init` option is only allowed for new runs. When the option is
used. A netCDF file, `INPUTDIR/REF_PBOT_FILE`, is written to store the
variable `REF_PBOT_VARNAME`, which can be used for restart runs.

* For restart runs, only `file` option is allowed. If `init` is
specified, a warning message is given and `SAL_REF_PBOT_CONFIG` is
changed into `file`.

* The check for restart runs does not work for unsplit mode, as restart
control structure is not accessible to unsplit dyncores. This will be
fixed later

* Note that any initial surface pressure from the atmosphere and
ice is not accounted for with this option.

  Revised the MOM_ice_shelf_diag_mediator code that creates the
MOM_IceShelf.available_diags file to follow the same pattern of logging as the
rest of MOM6.  The specific changes include moving the module name(s) to their
own line, adding static diagnostics to the variables that are logged in the
available_diags file (however, the ice shelf code does not have any static
diagnostics yet), logging the axes of the diagnostic, and adding cell metrics to
describe how a diagnostic might be reduced and logging these as well.  All
solutions are bitwise identical, but there are new entries and added information
in the MOM_IceShelf.available_diags files.

  Modified set_IS_axes_info() to follow the pattern from the rest of MOM6 and
take the units for the horizontal axes from the grid type.  In most cases, this
has the effect of changing the units from "degrees_E" and "degrees_N" to
"degrees_east" and "degrees_north" to follow CMIP protocols, and to mirror
changes that were made to the standard MOM6 framework code in 2017.  The names
of the gridspace axes were changed from `xB`, `yB`, `xT` and `yT` to `Iq`, `Jq`,
`ih` and `jh` to mirror the names in MOM6.  The long descriptions of these array
axis variables were also updated for clarity when `USE_INDEX_DIAGNOSTIC_AXES` is
true.  All solutions are bitwise identical, but there are changes to the
metadata and axis variable names in some output files, and some runtime
parameters that had been read from both `set_IS_axes_info()` and
`set_grid_metrics()` in MOM_grid_initialize.F90 are now only read in the latter
routine.  As a result of these changes, `set_IS_axes_info()` no longer takes a
`param_file_type` argument and the grid argument is now `intent(in)`.

  Added a missing dimensional scaling argument in the test for when to write
harmonic analysis diagnostics in HA_accum.  Also used the new unscale argument
in 2 calls to real_to_time in HA_init and used the new scale argument in one
call to time_to_real in HA_accum.  Answers and diagnostics are bitwise identical
when scaling is not applied, but this does correct a problem where some
diagnostics may have been written at the wrong time previously when dimensional
rescaling of time is being used.

  Edited the conversion factors for the 'ale_u2' and 'ale_v2' diagnostics to use
an equivalent form that makes it a little more obvious that the conversion
factors correspond to the declared units.  All answers are bitwise identical.

  This commit clarifies a number of comments in the MOM_diag_mediator.  It also
corrects multiple instances of spelling errors or non-standard use of white
space in comments.  There was also some rearrangement of module use statements
and the elimination of other unneeded module use calls.  Only comments or
spacing are changed, and all answers are bitwise identical.

A locally allocated array was uninitialized which is harmless for most configurations
that use time-filtering in the vertical grid generation. The time-filter uses land
masks and sets things to zero on land. Without the time-filter, we encountered non-
reproducible behavior on the root PE.

I also de-allocated a local variable (Rb) used on another if-branch in the same routine.

This commit refactors update_OBC_segment_data using the previously
introduced new structure segment%field component. Specifically, the code
that calculates normal_vel, normal_trans, normal_trans_bt,
tangential_vel gradient and external values of tracer / thickness
reservoir, is updated. The changes include,

* The previous loop + if structure is replaced by direct indexing of the
physical fields.
* Tidal velocity and SSH code is refactored to cji structure (c for
constituent), rather than an embedded if-branch within [ji] loop of a
[ij]c structure.
    * As a result, segment has three new 2D arrays, tidal_v[nt] and
    tidal_elev.
* Introduce orientation-agnostic looping indices [ij][se]_seg, currently
only used for tracer cell fields.
* Add tr_index in OBC_segment_data_type (segment%field) to facilitate
reservoir update.

* Add Newton iterations for the ice-shelf velocity solution

The previous Shallow Shelf Approximation solution used an iteration-on-viscosity scheme with Picard
iterations only. This commit adds the capability to change to Newton iterations once a certain
tolerance is met, which should reduce the number of iterations needed to reach convergence. This
tolerance to change to Newton is set by the new parameter NEWTON_AFTER_TOLERANCE. This parameter
defaults to the same value of ICE_NONLINEAR_TOLERANCE so that Newton is not called by default. The
Newton scheme is paired with the Eisenstat & Walker 1994 approach for adapting the CG tolerance to
help the Newton scheme achieve quadratic convergence without over-tightening at later Newton
steps. To use the adaptive CG tolerance, set NEWTON_ADAPT_CG_TOLERANCE=True (default is True).

* Newton iteration PR fix

Style changes, parentheses, removed string comparison from loops

* Pass_var optimization and FMAs

Made some changes to pass_var 'complete' arguments for Newton-related
variables. Also removed some Newton-related pass_vars and a
CS%doing_newton=.false. that were made before entering the the SSA
solution loop (at this point, Newton is already false, and there is no
need to update Newton-related halos yet). Added parentheses for
rotational symmetry with FMAs.

  Replaced "end if" with "endif" in 35 places in 10 files and "end do" with
"enddo" in 12 places in 4 files to comply with the MOM6 style guide.  In
addition, in 23 places in 7 files, code like `if(test)` were replaced with
`if (test)`, also to align with the MOM6 style guide.  All answers are
bitwise identical.

  Removed 41 trailing semicolons (MOM6 is in Fortran, not C).  In 406 other
lines of code the white space around semicolons was standardized to ' ; ' to
match the other ~14000 lines of code with semicolons joining Fortran statements
on the same line.  In a few cases there were also some improvements to indents
related to these changes.  With these changes, all semicolons that could be
replaced with a newline and still have correct Fortran follow the ' ; '
pattern.  These changes are scattered across 86 files.  Most of these changes
involve only white space, and all answers are bitwise identical.

Create a new subroutine read_OBC_segment_data out of a very chunky
update_OBC_segment_data. The new subroutine enables better control on
whether the update is associated with new external data through IO or
just recalculation (with evolved model internal state, i.e thickness,
tides).

Misc
* Use orientation-agnostic loop indices for segment%h and segment%htot
* Add [ij]_offset_in to loop over interior cells
* Field segment%Cg, which is never used, is now removed.

Create a new subroutine initialize_OBC_segment_reservoirs to initialize
thickness and tracer reservoir of OBC segments. This part of code was
original part of update_OBC_segment_data, but it is only called during
initialization.

This commit fixes a bug in particle advection of particles that are advected by uh and vh. When the tracer timestep doesn't match the dynamics timestep, these particles need to be advected using the tracer timestep. The previous version of the code sometimes leads to unphysical particle movement, and this commit fixes the issue.

  Added the new interfaces `symmetric_sum()` and `symmetric_sum_unit_tests()` to
the MOM_array_transform module, and added a call to `symmetric_sum_unit_tests()`
to `unit_tests()`.  The two functions wrapped by the `symmetric_sum()` interface
are intended to be a general template for rotationally symmetric sums, but they
demonstrably work as intended (as demonstrated by the passing unit tests) for
any rectangular arrays.  The new driver in test_MOM_array_transform.F90 tests
these new arrays.  Although `symmetric_sum()` is not used yet in the MOM6 code
apart from the tests, it was developed to address problems with failures of
rotational symmetry in the diagnostic downscaling routines, and will be used
there soon.  All answers are bitwise identical, but there are two new publicly
visible interfaces.

Resolves merge conflict due to new PLM_WLS test.

This commit refactors normal velocity/transport calculation in
update_OBC_segment_data, using the orientation-agnostic loop ranges.
Note that there is still an orientation-branch because the grid spacing
G%dyCu and G%dxCv are different.

* The field segment%h is removed, as we can just use h directly.
* Comments are added to better illustrate the subroutine structure

Remove inaccurate warning message in logs when OBC_TEMP_SALT_NEEDED_BUG
is false but temperature and salinity are needed and provided.

This commit fixes the bug that OBC segment data does not update if only
value mode is used for all segments. Previously, there were two issues,

1. Subroutine update_OBC_data, which is a wrapper over subroutine
update_OBC_segment_data is only called if global flag OBC%update_OBC is
True. This flag, however, is only turned on if there OBC is specified by
files.
2. Inside of update_OBC_data, subroutine update_OBC_segment_data is only
called if global flags OBC%any_needs_IO_for_data or
OBC%add_tide_constituents is True, which again prohibits OBC update with
"value".

Example: if normal velocity is specified by value, the transport is
calculate once during initialization using initial thickness. And the
transport is never updated through the run. And the model will restart
with a different transport.

A new runtime parameter OBC_VALUE_UPDATE_BUG is added to fix/recover
this bug. In addition, OBC%update_OBC is set to True if
initialize_segment_data is called. a separate procedure is added to set
any_needs_IO_for_data, for clarity.

  Added diag_mediator capabilities that had been developed and used in the SIS2
and MOM_ice_shelf versions of the diag mediator, in preparation for some of
these separate versions to be combined into a single version.  This includes
adding an overload to the register_scalar_field() interface so that scalars can
be registered with or without providing a null axis, retaining the existing
interface while adding the one that is widely used with the ice shelf code.  The
previous routine register_scalar_field() is now register_scalar_field_CS(), but
it has the same public interface.

  This commit also adds the new public interface MOM_diag_send_complete() that
calls the FMS routine that flushes the IO buffers.  Although this call can be
too expensive for production runs, it can be very useful for debugging.

  Internally, register_diag_field() now checks for standard axes for 2-d fields,
analogously to what was previously done for 3-d fields, thereby avoiding the
allocation of a separate axis group for each 2-d diagnostic.  This has been the
case for ice-shelf code for some time, and it will save memory will giving
identical results.

  The testing for incompatible mask and array sizes was consolidated into a
single call each in post_data_2d_low() and post_data_3d_low(), paving the way
for adding the option of masking static fields.

  All solutions are bitwise identical and diagnostic output (including metadata)
is identical to what was previously output, but there are two new public
interfaces.

  Renamed the q-point index space axis names to be uniformly Iq and Jq, rather
than having names that change depending on whether or not symmetric memory is
being used.  This commit also changes the long names of the grid space axes to
be consistent between the x and y axes and to mirror those that have been
adopted in the ice shelf code.  The order of the logical tests for index space
axes and symmetric memory around the calls to diag_axis_init were swapped to
that corresponding cases are next to each other, perhaps making any
inconsistencies more obvious.  There were also minor changes to the code setting
the axis label variables to follow the MOM6 case convention for indices and the
spacing conventions for assignments, as described in the MOM6 style guide.  All
solutions are bitwise identical, but there are changes to the names of the axis
label variables and their long descriptions when USE_INDEX_DIAGNOSTIC_AXES is
true.

  Corrected the cell_methods for 8 families of tracer content tendency
diagnostics.  Without this change, the 2-fold downscaled tendencies are 4 times
too large, and the 3-fold downscaled tendencies 9 times too large.  The
misleading comment describing the local `conv_units` string variable in
register_tracer_diagnostics() was also fixed.  All solutions and standard
diagnostics are bitwise identical, but there are changes in the values of
downscaled diagnostics and changes in metadata.

Makedep preprocessor parsing assumed define macros of the form `#define
macro(X) some-macro` but did not support empty macros of the form
`#define macro(X)`.

Preprocessor flag macros were assumed to be of the form `-DMACRO=VALUE`,
and did not support empty macros `-DMACRO`.

This patch addresses both issues.

…_vars

- To reproduce across restart files requires an additional halo update on diffu and diffv in MOM_dynamics_split_RK2.F90
- To reproduce across layouts with vertex shear option requires a halo update on thickness before remapping the viscosities on the vertex points.

- Setting BBL_USE_TIDAL_BG = True causes the model to fail dimensional consistency test.  If you set the spatial scale as m_to_L instead of m_to_Z in MOM_set_viscosity it fixes the issue.
- Fixing documentation fo L/T scaling of tideamp
  - Change suggested Hallberg

#1079)

This PR corrects an implementation of a finite-element solver for ice-sheet/shelf velocity. It computes a Jacobian weight J_q in shape functions for quadrature points and applies it to integrals over cells. This correction is significant when the shape of the cell is not rectangular. The same corrections are applied for subgrid parameterizations.

Additionally, it fixes missing deallocations of several arrays.

Changes made:
 - nitialize_ice_shelf_dyn: allocates a new variable `Jac`
 - bilinear_shape_fn_grid: computes a Jacobian determinant (`Jac`) at each quadrature point
 - CG_action: uses `jac_wt= Jac * IareaT` for computing the ice deformation and basal traction terms
 - matrix_diagonal: uses `jac_wt= Jac * IareaT` for computing diagonal terms
 - CG_action_subgrid_basal: computes `jac_sub_wt` using new arguments `dxCv_S`, `dxCv_N`,`dyCv_E` and `dyCv_W` that determine cell-edge spacing
 - ice_shelf_dyn_end: deallocates `Jac`, `CS%sx_shelf`, `CS%sy_shelf`, `CS%u_flux_bndry_val`, `CS%v_flux_bndry_val`,`CS%calv_mask`, `CS%Phi`, `CS%Phisub`, `CS%PhiC`.

Copilot (model Claude Sonet 4.6) assisted with this PR.

…1089)

* Vendor sphinxcontrib-autodoc_doxygen into docs/_ext/

This is piece 1 of the sphinx toolchain upgrade.  The docs build
previously depended on four forks which haven't been updated since 2020.

The chain pinned the entire build to Sphinx 3.2.1, requiring a growing
list of transitive version ceilings (jinja2<3.1, sphinxcontrib_*<1.0.x,
alabaster<0.7.14, setuptools<82).

This commit handles the sphinxcontrib-autodoc_doxygen dependency. The
upstream (rmcgibbo) has been dormant since June 2021; the fork was
effectively a rewrite (~90% of xmlutils.py changed, both documenters
replaced, the whole domain switched from cpp to f). Monkey-patching
the upstream package was rejected because there was nothing stable to
patch against. Instead the fork's code is vendored in-tree at
docs/_ext/autodoc_doxygen/ so it can be debugged, blamed, and edited
like any other project source.

Changes vs the fork at tag 0.7.13:

 - Renamed doxynamespace.rst template to doxymodule.rst and dropped
   the unused C++ doxyclass.rst (DoxygenClassDocumenter was already
   unregistered in the fork).
 - Removed ~80 lines of dead code: commented-out pdb.set_trace lines,
   the unregistered DoxygenClassDocumenter class, the dead
   visit_ref_angus alternative, _import_by_name_original, and the
   unused try-import of flint.
 - Dropped Python 2 compatibility cruft: from __future__ imports,
   from six import itervalues (replaced with dict.values()).
 - Ported four Sphinx 3 -> 8 API changes:
     * DoxygenAutosummary.get_items now passes self.bridge (the
       DocumenterBridge) to documenter constructors instead of self,
       so Documenter.__init__ finds directive.genopt.
     * self.warn / self.directive.warn calls (removed in Sphinx 8)
       replaced with sphinx.util.logging.getLogger(...).warning.
     * Wrapped env.doc2path in os.fspath to silence the
       RemovedInSphinx90Warning about str paths.

docs/conf.py: prepend _ext to sys.path and rename the extension in
the extensions list from sphinxcontrib.autodoc_doxygen to
autodoc_doxygen.

docs/conf.py also carries a monkey-patch for a parallel-build bug in
upstream VACUMM/sphinx-fortran's FortranDomain.merge_domaindata. The
upstream method references an undefined `outNames` (typo for
`ourNames`) and also unpacks the wrong tuple shape from modules and
objects dicts. With -j > 1 the merge silently fails and the f domain
ends up empty, losing every Fortran cross-reference. The patch is
heavily commented with a TODO to submit an upstream PR and drop the
workaround during piece 2 of the upgrade.

docs/requirements.txt: drop the sphinxcontrib-autodoc_doxygen git
dependency. The other three forks remain pending pieces 2, 3, 4.

Validation: built against a new venv.sph8 (Python 3.11 + Sphinx 8.2.3 +
sphinx-rtd-theme 3.1.0 + sphinxcontrib-bibtex 2.6.5 + upstream
VACUMM/sphinx-fortran master). `make html` produces 220 warnings vs
224 on the Sphinx 3 baseline. The f domain indexes 48 modules and
538 objects. Equation pages are effectively byte-identical (off only
by Sphinx 8's section-id hoist optimization and the pilcrow-as-CSS
theme change). Module pages are slightly larger because more
cross-references now resolve than in the baseline build.

* Swap sphinx-fortran fork for upstream VACUMM pinned commit

Piece 2 of the sphinx toolchain upgrade.

The jr3cermak/sphinx-fortran@1.2.2 fork existed to add Sphinx 3 API
compatibility in 2020. Upstream VACUMM/sphinx-fortran has continued
evolving since then (Sphinx logging API migration, parallel read
support, modern setuptools, removal of the future dependency) and now
works with Sphinx 8 directly. Upstream has not cut a PyPI release past
1.1.1 so we pin to a specific master commit for reproducibility.

The pinned commit ships a broken FortranDomain.merge_domaindata that
loses every f-domain object in parallel builds; the monkey-patch in
conf.py from the previous commit handles that. We will revisit that
patch and the pin together once an upstream PR is merged (see the
TODO(piece-2) marker in conf.py setup).

Validation: rebuilt docs with the pinned commit installed into
venv.sph8. `make html` produces 220 warnings (identical to the
pre-swap sph8 build), the f domain indexes 48 modules and 538 objects,
mom.html has 36 internal cross-reference links, and f-modindex.html is
generated. No behavior change from the sph8 build validated in the
previous commit.

* Drop sphinx fork for stock Sphinx 8

Piece 3 of the docs toolchain upgrade plan.
This commit removes the last need for the patched sphinx fork and
collapses the long list of transitive version ceilings it forced.

Changes:

 - requirements.txt: drop git+https://github.com/jr3cermak/sphinx.git
   @v3.2.1mom6.4. Pin sphinx>=8,<9 from PyPI. Drop the eight
   ceiling lines (jinja2<3.1, sphinxcontrib_*<1.0.x, alabaster<0.7.14,
   setuptools<82.0.0) and the workaround comment around sphinx-fortran's
   broken requirements.txt; those existed only because Sphinx 3.2.1 dragged
   them in. Add lxml as an explicit dep (it is required by the vendored
   _ext/autodoc_doxygen extension to parse Doxygen XML). Keep six because
   the pinned sphinx-fortran commit imports it at module load time.

 - conf.py: add a monkey-patch for sphinx.util.math.wrap_displaymath that
   replicates the functional changes in the jr3cermak fork. The patch
   detects parts containing `\begin{equation}`, `\begin{eqnarray}`, or
   `\begin{align}` and emits them verbatim (no outer
   `\begin{equation}\begin{split}...\end{split} \end{equation}` wrapping).
   Without this, MOM6's math-heavy sources produce nested LaTeX
   environments that pdflatex chokes on. The patch is heavily commented
   with the rationale and a TODO marker for a possible upstream
   contribution. Verified by unit-testing the function directly (plain
   math still wrapped, explicit-environment math passed through verbatim)
   and by inspecting the generated _build/latex/ MOM6.tex: 0 double-wrapped
   equation/eqnarray/align blocks, 100 + 18 + 17 explicit environments
   preserved as top-level blocks.

 - conf.py: tighten exclude_patterns. Sphinx walks the entire source
   tree by default, which previously caused it to descend into local
   virtualenvs (`venv.sph3/`, `venv.sph8/`, `venv-3.11/`, etc.) and
   pick up LICENSE.rst, README.rst, and stray autosummary template
   files from inside site-packages. Each of those generated a "document
   isn't included in any toctree" warning and a corresponding stray
   .html file in the build output. The new exclusions cover venv*,
   venv-*, _build.*, and _ext/*/templates. Final build: 90 HTML files
   (down from 148, the lost 58 were all venv junk and template
   scaffolding) and 122 warnings (down from 220, with the same project-
   internal warnings as before plus the venv noise removed).

 - Makefile: add the equation post-processing hook to the html target,
   gated on UPDATEHTMLEQS=Y. The jr3cermak/sphinx fork carried this
   as a patch to sphinx/cmd/build.py so it ran after every sphinx-build
   invocation. Now invoked from the Makefile so it works against
   stock upstream Sphinx. The existing nortd target already had the
   same hook with different arguments (-p APIs -b doxygen vs -p html
   -b sphinx).

The forked sphinx-fortran fix and the upstream PR for the math wrapping
change are deferred to a follow-up.

Validation: built from a fresh venv (venv.sph8.fresh) created by
`python3.11 -m venv` and `pip install -r requirements.txt`. The
install resolves cleanly with no manual intervention. `make html`
succeeds with 122 warnings, indexes 48 modules and 538 f-domain
objects, generates f-modindex.html, and produces 36 internal cross-
reference links on mom.html. `make html UPDATEHTMLEQS=Y` runs the
post-processing script successfully. `make latex` produces 0
double-wrapped math environments.

* Drop flint dependency from docs build

Piece 4 of the docs toolchain upgrade plan.
With this commit the dependencies in docs/requirements.txt are only on
maintained upstream sources.

flint (marshallward/flint, vendored as jr3cermak/flint at 0.0.1) was
historically pulled into the docs build with the intent of patching
up Doxygen's incomplete parsing of Fortran functions with `result()`
clauses. The fork's sphinxcontrib-autodoc_doxygen imported it at
module load time but never actually called it; the import was a
no-op behind a `try: import flint; except: pass`. We removed that
import as part of piece 1's dead-code cleanup, so the vendored
extension no longer touches flint at all.

Verified empirically: uninstalled flint from venv.sph8.fresh and
rebuilt. The build is byte-identical to the previous one with flint
installed (122 warnings, 48 f-modules, 538 f-objects, 36 internal
cross-reference links on mom.html, f-modindex.html present, 90 html
files). Nothing in the current pipeline depends on flint.

The two surviving "flint" mentions in docs/_ext/autodoc_doxygen/ are
both in comments describing the original intent, not live code; they
are left in place as historical context. If a future enhancement
wants to actually fix Doxygen's result() clause parsing, that should
be a separate piece of work and would not necessarily resurrect
flint specifically.

docs/README.md is updated to drop the flint bullet from the
requirements list and the flint row from the credits table. Other
stale entries in that section (sphinxcontrib_autodox-doxygen, the
sphinx and sphinx-fortran fork rows) remain and will need a broader
README sweep in a follow-up.

* Sphinx upgrade follow-up fixes

Consolidates four small fixes that followed the Sphinx 8 upgrade:

- Switch all Makefile targets from `sphinx-build -b` to `-M` for
  consistent parallel build support.
- Update docs/README.md to match the modernized toolchain (remove
  references to the four sphinx* forks, update requirements and
  credits sections).
- Fix crash in autodoc_doxygen visit_image on empty <image>
  elements (node.text is None when doxygen produces <image/>
  without caption text).
- Add __pycache__ and *.pyc to docs/.gitignore so the vendored
  extension's bytecode does not appear in git status.

* RTD build infrastructure

Consolidates four fixes for the Read the Docs build environment:

- Override RTD's default sphinx runner with a build.jobs.build.html
  entry in .readthedocs.yml that runs sphinx-build -M html -j auto,
  since RTD's high-level sphinx: key has no parallelism option.
- Make doxygen_xml path resolution robust to cwd changes: resolve
  relative paths against app.confdir rather than the ambient cwd,
  which differs between local builds and RTD.
- Switch Makefile html target from -j 4 to -j auto to match RTD.
- Enable html_static_path = ['_static'] in conf.py so custom CSS
  files (autodoxysource.css) are copied into the build output.
  Previously commented out, which meant app.add_css_file() added
  the <link> tag but the file was never deployed.

* Fix O(N^2) XPath scans in autodoc_doxygen

Two instances of the same class of bug in the vendored
autodoc_doxygen extension, both scanning the entire merged doxygen
XML tree (1230 files, 109 MB with programlisting) on every call:

1. scanNode used `node.xpath('//latexonly')` etc. In XPath, `//`
   starts from the document root, not from `node`. Since the
   extension merges all XML into one tree, every call scanned
   the whole tree. Fix: `//` -> `.//` (descendant-of-self).
   Recovery: 6m33s -> 48.8s wall clock at full input scale.

2. visit_ref used `get_doxygen_root().findall('.//*[@id=X]')` to
   resolve each prose <ref> by linearly scanning every element in
   the merged tree. With XML_PROGRAMLISTING=YES the tree tripled
   in size, making this the dominant cost (250s of 911s serial).
   Fix: lazy {id: element} dict built once on first use, O(1)
   lookup thereafter. Recovery: 911s -> 676s serial.

* Add source code browser to Sphinx docs

Generate one HTML page per Fortran source file with syntax
highlighting and clickable cross-references. Clicking an
identifier in the source jumps to its API documentation page;
each API entry gets a [source] link back to the source listing.

Implementation:
- Enable XML_PROGRAMLISTING in Doxyfile_rtd (standalone commit
  originally, now folded in).
- New autodoxysource directive in _ext/autodoc_doxygen/ that walks
  doxygen's <programlisting> XML, emitting per-line anchors,
  CSS-classed highlight spans, and pending_xref nodes for
  identifiers.
- Stub generator produces 329 :orphan: source pages under
  api/generated/source/.
- [source] links added to DoxygenMethodDocumenter,
  DoxygenTypeDocumenter, and DoxygenModuleDocumenter.
- CSS styled to match the sphinx_rtd_theme's code blocks (font
  stack, size, colors, line gutter).
- Node count optimization: coalesce text runs in _walk_highlight
  (10-30 nodes/line -> 3-5), set support_smartquotes=False on the
  source container to skip the smartquotes transform. Combined
  with the visit_ref fix, total build time with the source browser
  is 115s parallel (vs 360s before optimization).

* Improve API documentation

Consolidates five improvements to the rendered API pages:

- Show Fortran type declarations for function/subroutine
  parameters and derived-type members. Types are rendered as
  inline code (e.g. ``real, intent(in), optional``) by looking
  up the <param> elements on the parent <memberdef>.
- Move [source] links from the top of each entry to the bottom,
  so the description and parameters appear first.
- Fix functions (as opposed to subroutines) not being registered
  in the sphinx-fortran domain. format_name() was prepending the
  return type (e.g. "real") which sphinx-fortran's f_sig_re regex
  cannot parse, leaving function entries without anchors.
- Add Functions and Source Files index pages to the API Reference,
  with cross-reference links to module pages and source browser
  pages respectively.
- Exclude upstream CVMix sources from doxygen input to avoid
  warnings from undocumented third-party code.

* Fix multi-line math blocks losing indentation in LaTeX

Multi-line math content (e.g. \begin{align}...\end{align}) from
doxygen XML was emitted as:

    .. math:: \begin{align}
    \mathbf{v}
      = \mathbf{u} + ...
     \end{align}

RST requires directive body content to be indented. The
continuation lines at column 0 were parsed as regular text by
docutils, producing garbled LaTeX output with "Runaway argument"
and "Missing $" errors — 164 LaTeX errors on the Notation page
alone.

Fix: in visit_formula, detect multi-line math text and emit it
as a properly indented block:

    .. math::

       \begin{align}
       \mathbf{v}
         = \mathbf{u} + ...
       \end{align}

Single-line math is unchanged (stays on the .. math:: line).

After fix: `make latexpdf` LaTeX errors drop from 164 to 10.
The remaining 10 are unrelated (9 Unicode Greek characters in
source comments that pdflatex cannot render, 1 duplicate label).

Under optimization -O2 with ifort we had some floating-point divergences
in the new Recon1d schemes. Two schemes diverged from their counterparts
that use the old "select" pathways. One scheme failed self-consistency
tests due to inlining allowing some aggressive optimization to occur.

Recon1d_PPM_H4_2019:
- Add explicit parentheses in end_value_h4() to enforce evaluation order
  and ensure bit-reproducibility across compilers.

Recon1d_PPM_hybgen:
- Refactor reconstruct() into a cleaner two-pass approach using new
  private helpers (hybgen_ppm_coefs, bound_edge_values,
  check_discontinuous_edge_values) copied from phased-out modules to
  avoid external dependencies. Bit-for-bit with the old PPM_HYBGEN scheme.
- Add average() override that clamps sub-cell averages to
  [min(ul,ur), max(ul,ur)], reproducing force_bounds_in_subcell behaviour.

Recon1d_PLM_WLS:
- Fix LS_error(): remove 0.5/0.25 factors and reformulate as squared
  deviation from the optimal slope.
- Also updated Doxygen math which had some errors.
- In check_reconstruction(), copy state via assignment (perturbed = this)
  since an inline call to %reconstruct() on same data was yielding LSB
  differences.

* More fixes to achieve consistent rotation with FMAs

* New Krylov methods for ice-shelf velocity

Added MINRES and CR methods to solve the shallow shelf approximation.
Further optimized CG. The new CG is answer changing, but only because
IDIAGU/IDIAGv reciprocals are now procomputed to replace per-iteration
divisions with multiplications. MINRES and CR can be faster than CG.
Also added CG_HALO_SHRINK parameter to disable the halo-shrinking
scheme for CG; combined with the other changes, this reduces the
number of pass_vector calls to once per iteration, rather than 4
pass_vector calls every nihalo iterations. Disabling the
halo-shrinking can be faster for certain combinations of grid
properties (halo size, symmetric/nonsymmetric, grid size).

* Remove unneeded line

Removed a line that caused an error when compiling in debug mode,
where a variable was being set to the value of another variable that
was not set yet. This line was not needed.

* Add wrapper for ice-shelf Krylov methods

Added a wrapper subroutine that performs the shared setup for each of
the Krylov methods, calls the selected Krylov method, and assigns
BCs. Also: made sure the number of iters is returned from each method
in the unlikely case of an early exit; slight change to where CG
convergence is tested to avoid evaluating extra unneeded code.

* Fix ice-shelf inner solver units. Change exit point for MINRES for efficiency.

* Fix units for residual scaling within inner solvers

* Fix formatting